home *** CD-ROM | disk | FTP | other *** search
/ Amiga Plus 2000 #4 / Amiga Plus CD - 2000 - No. 4.iso / PowerPC / Dev / PPCRelease / Libmfd / s_asinh.c < prev    next >
Encoding:
C/C++ Source or Header  |  1997-01-02  |  1.6 KB  |  62 lines
  1.  
  2. /* @(#)s_asinh.c 1.3 95/01/18 */
  3. /*
  4.  * ====================================================
  5.  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  6.  *
  7.  * Developed at SunSoft, a Sun Microsystems, Inc. business.
  8.  * Permission to use, copy, modify, and distribute this
  9.  * software is freely granted, provided that this notice 
  10.  * is preserved.
  11.  * ====================================================
  12.  */
  13.  
  14. /* asinh(x)
  15.  * Method :
  16.  *    Based on 
  17.  *        asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
  18.  *    we have
  19.  *    asinh(x) := x  if  1+x*x=1,
  20.  *         := sign(x)*(log(x)+ln2)) for large |x|, else
  21.  *         := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
  22.  *         := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))  
  23.  */
  24.  
  25. #include "fdlibm.h"
  26.  
  27. #ifdef __STDC__
  28. static const double 
  29. #else
  30. static double 
  31. #endif
  32. one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
  33. ln2 =  6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
  34. huge=  1.00000000000000000000e+300; 
  35.  
  36. #ifdef __STDC__
  37.     double asinh(double x)
  38. #else
  39.     double asinh(x)
  40.     double x;
  41. #endif
  42. {    
  43.     double t,w;
  44.     int hx,ix;
  45.     hx = __HI(x);
  46.     ix = hx&0x7fffffff;
  47.     if(ix>=0x7ff00000) return x+x;    /* x is inf or NaN */
  48.     if(ix< 0x3e300000) {    /* |x|<2**-28 */
  49.         if(huge+x>one) return x;    /* return x inexact except 0 */
  50.     } 
  51.     if(ix>0x41b00000) {    /* |x| > 2**28 */
  52.         w = __ieee754_log(fabs(x))+ln2;
  53.     } else if (ix>0x40000000) {    /* 2**28 > |x| > 2.0 */
  54.         t = fabs(x);
  55.         w = __ieee754_log(2.0*t+one/(sqrt(x*x+one)+t));
  56.     } else {        /* 2.0 > |x| > 2**-28 */
  57.         t = x*x;
  58.         w =log1p(fabs(x)+t/(one+sqrt(one+t)));
  59.     }
  60.     if(hx>0) return w; else return -w;
  61. }
  62.